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Abstract. We consider the head-on collision of equal-mass neutron stars boosted towards 
each other and we study the behavior of such systems near the threshold of black-hole 
formation. In particular, we confirm the previous findings by Q that a type-I critical 
phenomenon can be observed by fine-tuning the initial mass of the two neutron stars. 
f*^ , At the same time, we argue against the interpretation that the critical solution is not a 

*vj ■ perturbed spherical star and show instead that the metastable star corresponds to a (perturbed) 

equilibrium solution on the unstable branch of the equilibrium configurations. As a result, the 

^ ' head-on collision of two neutron stars near the critical threshold can be seen as a transition 

Q ' in the space of configurations from an initial stable solution over to a critical metastable one 

^~^ , which can either migrate to a stable solution or collapse to a black hole. The critical exponent 

~ ' , for this process shows a fine structure which was already observed in the case of the critical 

£*\J ■ collapse of scalar fields but never before for perfect fluids. 

O | PACS numbers: 04.25.Dm, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60Jd 
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1. Motivation and introduction 



In 1993, Choptuik (2| considered one-parameter families of solutions, S[P], of the Einstein- 
Klein-Gordon equations for a massless scalar field in spherical symmetry, such that for every 
J^ P > P*, S[P] contains a black hole and for every P < P*, S[P] is a solution not containing 

qv ' singularities. He studied numerically the behavior of S[P] as P — > P* and found that the 

^. , critical solution, <S[P*], is universal, in the sense that it is approached by all nearly-critical 

CS| ' solutions regardless of the particular family of initial data considered. He also found that S [P] 

r^. , exhibit discrete self-similarity and that, for supercritical solutions (P > P*), the mass of the 

black hole satisfies .A/bh = c\P — P*] 1 , with 7 being an universal constant, i.e. not depending 
^D . on the particular family of initial data. 

After Choptuik's seminal work, similar transitions were discovered for a wide range 
of systems, including massive scalar fields and ultra-relativistic fluids, see Q for a recent 
review. All these phenomena have the common property that, as P approaches P*, S[P] 
approaches a universal solution 5[P*] and that all the physical quantities of S[P] depend 
only on \P — P*\. In analogy with critical phase transitions in statistical mechanics, these 
transitions in gravitational collapse were later classified as "type-I" critical phenomena, with 
static or periodic critical solutions and discontinuous transitions in the vicinity of the critical 
point, or "type-II" critical phenomena, with self-similar critical solutions and continuous 
transitions in the vicinity of the critical solution |3). 

The study of critical phenomena in neutron-star (NS) collapse started with the work 
by H on radiation fluids and was later extended to more general ultra-relativistic equations of 
state (EOS) 15161 and ideal-gas EOS (2r[9). In all these studies the collapse was triggered using 
strong perturbations and a type-II critical phenomena was found. Type-I critical phenomena 
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in the collapse of unstable configuration under very small perturbations was instead studied 
only very recently II 1 Oil 1 1 H . 

The first study of critical phenomena the head-on collision of NSs was instead carried 
out by Jin and Suen in 2007 jTJ. In particular, they considered a series of families of equal 
mass NSs, modeled with an ideal-gas EOS, boosted towards each other and varied the mass 
of the stars, their separation, velocity and the polytropic index in the EOS. In this way they 
could observe a critical phenomenon of type I near the threshold of black-hole formation, 
with the putative solution being a nonlinearly oscillating star. In a successive work JT2), they 
performed similar simulations but considering the head-on collision of Gaussian distributions 
of matter. Also in this case they found the appearance of type-I critical behaviour, but also 
performed a perturbative analysis of the initial distributions of matter and of the merged 
object. Because of the considerable difference found in the eigenfrequencies in the two cases, 
they concluded that the critical solution does not represent a system near equilibrium and in 
particular not a perturbed Tolmann-Oppenheimer-Volkoff (TOV) solution QJ. 

In this paper we study the dynamics of the head-on collision of two equal-mass NSs using 
a setup which is as similar as possible to the one considered in JT|. While we confirm that 
the merged object exhibits a type-I critical behaviour, we also argue against the conclusion 
that the critical solution cannot be described in terms of equilibrium solution. Indeed, we 
show that, in analogy with what found in ifTTI . the critical solution is effectively a perturbed 
unstable solution of the TOV equations. Our analysis also considers fine-structure of the 
scaling relation of type-I critical phenomena and we show that it exhibits oscillations in a 
similar way to the one studied in the context of scalar-field critical collapse 11311141 . 

The remainder of this paper is organized as follows. In section [2] we describe the 
numerical settings of the simulations, the properties of the used initial data and we define 
some of the quantities that we analysed. In section [3] we show in details our results, while 
section 2] is dedicated to conclusions and discussion. Finally, presented in Appendix A is a 
study of the compactness of the metastable solution and a comparison with its "hoop radius". 
The work reported is closely related to the analysis carried in a companion paper (hereafter 
paper I) about critical phenomena in linearly unstable nonrotating NS models irTTI . Because of 
the logical affinity between the two works we will often refer the reader to the results obtained 
in ifTTTl . We use a spacelike signature (— , +, +, +), with Greek indices running from to 3, 
Latin indices from 1 to 3 and the standard convention for the summation over repeated indices. 
Unless explicitly stated, all the quantities are expressed in the system of dimensionless units 
in which c = G = M© = 1 . 

2. Numerical setup and initial models 

In what follows we briefly describe the numerical setup used in the simulations and the 
procedure followed in the construction of the initial data. 

2.1. Numerical Setup 

The numerical setup used in our simulations is identical to the one presented in IfTTI . In 
particular, we use the axisymmetric code Whisky2D, which has been described in detail 
in 1(131 (but see also the more extended discussion in paper I) and is based on the 3-dimensional 
code Whisky irT6ffl"8l . to solve in 2 spatial dimensions the full set of Einstein and relativistic- 
hydrodynamics equations. The axisymmetry of the spacetime is imposed exploiting the 
"cartoon" technique introduced in 1191 . while the hydrodynamics equations (which are written 
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in first-order form and thus do not require double spatial derivatives) are written explicitly in 
cylindrical coordinates (see ifTSll to read about the advantages of this formulation). 

All the simulations presented here have been performed using an ideal-gas EOS, p = 
(r — l)pe, where p is the baryonic density and e the specific internal energy in the rest frame 
of the fluiqj. The discretization of the spacetime evolution equations is done using fourth- 
order finite-differencing schemes, while we high-resolution shock-capturing (HRSC) methods 
with a piecewise-parabolic method (PPM) reconstruction for the hydrodynamics equations. 
The time-stepping is done with a third-order total-variation diminishing (TVD) Runge-Kutta 
scheme and despite the use of higher-order methods, the convergence order drops to about 
1 after the two NSs have merged and large shocks develop (see l20l for a discussion on the 
convergence order in relativistic -hydrodynamics simulations). 

The spatial discretization is done a grid with uniform resolution, which we have taken 
to be either h = 0.1 M or h = 0.08 M . Furthermore, as long as nonspinning NSs are 
considered, the head-on collision also possesses a symmetry across the plane midway between 
the two stars and orthogonal to the colliding direction. As a result, the problem needs to be 
solved only for one star and suitable boundary conditions be applied across the symmetry 
plane. The outer boundary of the computational domain is set at 60 M and thus rather 
close to the two stars. However, this is adequate since we are not interested here in extracting 
gravitational waveforms and since we have verified that the violation of the constraints at the 
outer boundaries are not larger than elsewhere in the computational domain. 

As a final remark we note that we use the same gauges as those employed in 12TI and thus 
the slicing is sufficiently "singularity-avoiding" that it is not necessary to perform an excision 
of the field variables when following the evolution of a supercritical solution. However, 
because of the very high-resolution used, the rapid growth of the rest-mass density is not 
compensated by the intrinsic numerical dissipation as instead happens in l2"TI or in Ifl8l . As 
result we excise the solution of the hydrodynamical quantities only as discussed in 11711221 in 
order to obtain a stable, long-term solution. 

2.2. Initial data 

Our initial data consists of two equal-mass nonrotating NSs having initial coordinate 
separation, computed as the coordinate distance between the two stellar centres, of 20 M . 
Following A], we construct these stars using a polytropic EOS, p = Kp T , with adiabatic 
exponent T = 2 and polytropic constant K = 80, which is equivalent to K = 0.00298 c 2 /p n , 
and where p n = 2.3 x 10 14 g/cm is approximately the nuclear density. The stars in the 
neighbourhood of the critical solution have a radius of R = 14.7 km, a baryon/gravitational 
(ADM) mass of 0.760/0.732 M . The main properties of the critical solution are reported in 
tableQ] where they are indicated as model "A". 

Note that the maximum baryon/gravitational mass for the chosen value of the polytropic 
constant is 1.609/1.464 M and thus the object produced by the collision will have a 
baryon/gravitational mass which is above such maximum mass. However, as we will discuss 
below and in contrast with the claim made in 1 1 1, this can still lead to equilibrium solution for 
a TOV star. 

The stars constructed in this way are then boosted towards each other along the z- 
axis with a velocity vq = 0.15, which is similar to the free-fall velocity as computed from 
the Newtonian expression for a point-particle. With these choices the only remaining free 
parameter needed to characterize the initial data is the central rest-mass density of the two 

| We recall that an ideal-gas EOS and a polytropic one are identical for fluids undergoing isoentropic evolutions. 
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Table 1, Equilibrium models used for the collision (model A) or that are discussed in figure|5] 
(models B and C). Listed are: the value of the gravitational (ADM) mass, the total rest-mass, 
the radius of the star, the central rest-mass density and the polytropic constant K. For each 
model we also report the maximum allowed ADM and rest-masses for a TOV having the same 
polytropic constant. 



Model 


Madm 


M b 


ft (km) 


Pc 


K 


A/ADM,max 


-"-*£>, max 


A 


0.732 


0.760 


14.761 


0.00058 


80 


1.464 


1.609 


B 


1.505 


1.514 


9.135 


0.00963 


145 


1.972 


2.166 


C 


1.460 


1.547 


19.003 


0.00055 


155 


2.038 


2.240 



NSs, p c , which we use as our critical parameter (We note that the value chosen for the initial 
velocity does influence qualitatively the results obtained and indeed it can act as a critical 
parameter in a sequence at constant p c Q]). 

2.3. Entropy of the critical solution 

A convenient way to study the properties of the critical solution is to characterize its 
thermodynamical properties and in particular its entropy. The basic idea is that the 
metastable (critical) solution is simply an equilibrium solution on the unstable branch of 
TOV configurations and thus that by measuring its entropy it is possible to relate it to the 
corresponding equilibrium polytropic model. Following l23l . we express the second law of 
thermodynamics as 

S = S + C v m(T) - R ln(p), (1) 

where So is an integration constant that we set to zero and where 

c " - (£). 

is the specific heat capacity at constant pressure, and R the gas constant. In the case of an 
ideal gas R is related to the specific heat capacities through the relation 

where h = 1 + e + p/p is the specific enthalpy. Recalling that e = CyT and T = C p /Cy, 
the second law of thermodynamics ([TJ for a polytropic EOS is simply given by 

S = c - ln (cm^T)) ' (4) 

or equivalently 

K = R cxp \Jj-} . (5) 

Stated differently, the polytropic constant reflects all of the changes in the entropy of the 
system, so that expression ([5} allows a simple connection between the entropy of the 
critical solution, which we measure as proportional to K = p/p r , and the properties of a 
corresponding equilibrium TOV model. 



Critical Phenomena in Neutron Stars II: Head-on Collisions 5 

In practice the polytropic constant can change enormously across the star especially after 
the collision and since we are interested only in global quantities we use a volume-averaged 
polytropic constant 

, x Lp/p r dV 

and perform the volume integration not across the whole star but over a volume Q, where 
the rest-mass density is larger than 10% of the initial central one. This choice removes the 
difficulties with possible divergences near the stellar surface and we have verified that is robust 
against different values of the threshold density. 

3. Results 

We next discuss the dynamics of the collision and the properties of the critical solution as they 
appear from our numerical simulations. 

3.1. Dynamics of the collision 

The basic dynamics of the process is rather simple to describe. As the two NSs are accelerated 
towards each other by the initial boost and the mutual gravitational attraction, they collide, 
leading to a merged object which is wildly oscillating and with a mass which is above the 
maximum mass of the initial configuration. Depending on whether the initial central density 
is larger or smaller than the critical one, the metastable solution either collapses to a black hole 
(supercritical solutions) or it expands to a new stable stellar solution (subcritical solutions). 

Before entering the details of the discussion it may be useful to remark that the simulation 
of the head-on collision of two NSs in the neighborhood of the critical solution is a very 
demanding calculation, even for modern general-relativistic hydrodynamical codes. This is 
because adaptive mesh refinements do not provide any significant speed up and, at the same 
time, rather high level of resolutions are needed to capture the dynamics faithfully. In spite of 
these computational difficulties mentioned above, using the Whisky2D code we are able to 
identify the critical value for p c with an accuracy of 11 significant digits, a level of precision 
never achieved before in the study of NS head-on collisions. More specifically, we have 
measured the critical central density 

p* = (5.790998966725 ± 0.00000000003) x 10~ 4 . (7) 

as the midpoint between the largest central rest-mass density among the subcritical models and 
the smallest central rest-mass density among the supercritical ones. In other words, binaries 
with p c > p* have been simulated to collapse to a black hole, while binaries with p c < p* c 
have been computed to expand to a stable star. 

As pointed out by Q], this value will ultimately depend on the numerical resolution used 
and the other discretization parameters. Nevertheless, given a set of initial data and numerical 
setup, O gives a precise measure of how close we are able to get to the critical solution. 

Compressed in figure Q] is a considerable amount of information about the criticality of 
the head-on collision. More specifically, we show the time evolution of the central rest-mass 
density of the heaviest subcritical model (black solid line) and lightest supercritical model (red 
dashed line) computed. Overall, we can distinguish three different phases of the dynamics. In 
the first one (marked with "A" in the figure) the central density increases from its initial value 
to a maximum one reached, which is attained when the stellar cores merge. 
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Figure 1. Evolution of the maximum rest-mass density of the most massive subcritical model 
with initial value p c = 0.000579099896670 (black solid line) the least massive supercritical 
model with initial value p c = 0.000579099896675 (red dashed line). Also highlighted are 
the four different phases of the dynamics. The first one corresponds to the initial configuration 
of the system (labelled as "A"). After the collision, a new metastable solution is created 
during which the central density exhibits violent oscillations (labelled as "B"); the subrcritical 
and supercritical solutions are essentially indistinguishable during this stage. Finally the 
subcritical and supercritical solutions separate, with the first one relaxing to a stable expanded 
configuration (labelled as "C"), while the second collapses to a black hole (labelled as "D"). 



The second phase (marked with "B") starts from t ~ 0.4 ms and is characterized by 
strong oscillations around the metastable equilibrium. During these first two phases the 
subcritical solution d the supercritical one are essentially indistinguishable, but at the end 
of the second phase, i.e. at t ~ 1.2 ms, the different nature of the two solutions emerges 
and the evolutions of the rest-mass density differ. More specifically, during the third phase 
(marked with "C" and "D"), the supercritical solution shows an exponential increase of the 
central rest-mass density as a result of the collapse to a black hole. The subcritical solution, 
on the other hand, shows a violent expansion and the central density settles to a value which 
is about one tenth of the maximum one attained during the second phase. Clearly, the most 
interesting part is obviously the one corresponding to the second phase, during which the 
merged object is a metastable solution in which the central density has strong, non-harmonic 
oscillations (see inset of figure[T]i. 

In order to better describe the dynamics of the system, we show in figures [2] and [3] 
the evolution of colour-coded countours of the rest-mass density of the supercritical and 
subcritical solutions, respectively, when shown at representative times in the (x, z) plane. The 
first row of panels in figure[2]shows the initial configuration of the system at time t = and a 
subsequent stage, at time t = 0.394 ms, corresponding to when the first maximum in the rest- 
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Figure 2. Isodensity contours in the (x, z) plane of the least massive supercritical solution 
(p c = 0.000579099896675). The corresponding times are shown at of each panel, while the 
color-code for the rest-mass density is indicated to the right. The isolines are shown for the 
values of p= 10 12 ' 6 , 10 12 ' 7 , 10 13 ' 2 , 10 13 ' 7 , 10 14 ' 2 and 10 14 ' 7 g/cm 3 . The second frame 
(t = 0.394 ms) is taken during the merge process. The next five frames illustrate the star 
during the metastable equilibrium. Finally the last frame (i = 2.759 ms) shows the solution 
after the formation of a black hole, whose apparent horizon is shown with a thick black dashed 
line. 
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Figure 3. The same as figure [2] but for the most massive subcritical model (p c = 
0.000579099896670). The first five frames are similar (although not identical) to the 
corresponding ones in figure [2] since they refer to the metastable evolution when the two 
solutions are essentially indistinguishable. The last frame (i = 2.759 ms) shows a new NS 
produced from the migration of the metastable object to the stable branch of the equilibrium 
solutions. 
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mass density is reached (cf. figureQ]). This time also represents the one at which the two stellar 
cores enter in contact and thus marks the beginning of the metastable equilibrium. During this 
stage, two strong shock waves propagate along the z-direction, ejecting part of the matter as 
shown in the third panel at time t = 0.788 ms. Most of this matter is still gravitationally 
bound and falls back onto the central object creating a new shock wave (cf. fourth panel at 
t = 1.182 ms). This process is then repeated multiple times and results in a sequence of 
bounces until the object finally collapses to a black hole, as shown in the last panel at time 
t = 2.759 ms and which marks the fate of the supercritical solution. 

Similarly, figure [3] reports representative stages of the evolution of the subcritical 
solution. A rapid inspection of the first five panels of figure|2]indicates they are very similar to 
the corresponding ones in figure [2] Indeed, the supercritical and subcritical solutions are the 
same to the precision at which we measure the critical solution [cf. (0] and we have reported 
the panels here exactly to remark the similarity during the first two stages of the evolution. 
However, being it a subcritical solution, the metastable evolution does not end with a black 
hole formation but, rather, with a new stable stellar solution. This is shown in the sixth panel 
of figure|3]and refers to a time t = 2.759 ms, after the metastable star has expanded violently 
and when it has reached a new quasi-spherical configuration. 

With the exception that we are able to get closer to the critical solution, much of what 
reported here confirms what found by Jin and Suen in |1|. In the following section, however, 
we discuss how and why our conclusions about the properties of the critical solutions differ 
from those discussed in [fl~) and subsequently in |[T2l . 

3.2. Nearly-critical solutions 

The scope of this section is to show that, contrary to what suggested in lfTl fT2l . the metastable 
object can be interpreted rather simply as the perturbation of a new equilibrium configurations 
of linearly unstable spherical stars. To provide evidence that this is the case, we have 
computed the evolution of the average entropy of the system for a subcritical solution as 
computed via the volume-averaged polytropic constant ©. This is shown in figure [4] where 
we report both (K) and the central density p c during the metastable equilibrium or stage "B" 
(left panel) and when the metastable solution has expanded to recover a stable solution or 
stage "3" (right panel). 

Clearly, the two panels show that two quantities are correlated and indeed in phase 
opposition - entropy increases when the density decreases and viceversa - as one would 
expect from the first law of thermodynamics 

TdS = pdV + dQ= pdV , (8) 

where the second equality comes from assuming adiabatic transformations. 

Using the results shown in figure |4] it is then possible to compute a time-averaged value 
of (K) and of the central density, i.e. (K), ~pZ, and thus equilibrium polytropic models with 
such polytropic constant and central density. These equilibrium models can be constructed 
either relative to the metastable stage or relative to the final stable stage of the subcritical 
solutions. 

The results of this procedure are summarized in figure [5] which reports the equilibrium 
sequences of TOV stars with fixed polytropic constant in a standard (p c , Madm) plane. In 
particular, the black solid line refers to a sequence with K = 80 and the black triangle 
therefore marks the initial critical solution (i.e. model "A" in table [TJ. Similarly, the red 
long-dashed line refers to a sequence with a polytropic constant (K) = 145, which therefore 
coincides with the time averaged value of K during the metastable state and as deduced from 
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Figure 4. Left panel: Evolution of the central rest-mass density and of the effective polytropic 
constant (K) for the most massive subcritical model during the metastable equilibrium phase 
of the dynamics. Right panel: The same as in the left panel but during the relaxation to a stable 
configuration. 



figure [4] The red square shows therefore the equilibrium model having as central density the 
time-averaged central density of the metastable solution {i.e. model "B" in table [TJ. Finally, 
the blue dashed line refers to a sequence with a polytropic constant (K) = 155, which 
therefore coincides with the time averaged value of A" during the stable stage of the subcritical 
solutions. (c/figure|4|i. The blue pentagon shows therefore the equilibrium model having as 
central density the time-averaged central density of the stable solution {i.e. model "C" in 
table [1). It is important to remark that because models "A, B, C" are determined after fixing 
the polytropic constant and the central rest-mass density, they are not guaranteed to have 
the same total baryon mass. In practice, however, they do have the same rest-mass with a 
precision of 0.4% (of course 2Mb,A — Mb,B — Mb,c)- This is not a coincidence but a clear 
evidence of the common link among the three models. 

What is eloquently shown in figure [5] can also be stated summarized as follows. The 
head-on collision of two NSs near the critical threshold can be seen as series of transitions in 
the space of configurations from an initial stable model "A" over to a metastable model "B" 
which has the same rest-mass but larger gravitational mass as a result of the conversion of 
the kinetic energy into internal energy via large shocks. Because model "B" is on the linearly 
unstable branch of the equilibrium configurations, it can exhibit a critical behaviour (this was 
shown in great detail in paper I) and thus subcritical solutions will expand and move the stable 
branch of equilibrium solutions (model "C"), while supercritical solutions will collapse to a 
black hole (solution "D"). 

In the light of this interpretation, the conclusion drawn by Jin and Suen [TJ that the 
merged object is far from being a TOV because it promptly collapses even though its total 
rest-mass is smaller than the corresponding maximum mass, does not appear to be the correct 
one. Indeed, the transition highlighted in figure [5] clearly shows that even a sub-massive TOV 
can be brought over the stability threshold to collapse to black hole by simply increasing its 
gravitational mass, namely by increasing its internal energy. 

Additional evidence that the merged object is indeed a perturbed TOV comes from 
analyzing the oscillation frequencies measured over the metastable stage. Despite the latter 
is rather short and the eigenfrequencies are consequently not very accurate, they agree well 
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Figure 5. Equilibrium sequences of TOV stars with fixed polytropic constant in a standard 
(PcMadm) plane. The black solid line refers to a sequence with K = 80 and the black 
triangle marks the initial critical solution (i.e. model "A" in tablefT}. The red long-dashed line 
to a sequence with (K) = f45 and the red square shows the equilibrium model having as 
central density the time-averaged central density of the metastable solution (i.e. model "B" 
in table[TJ. Finally, the blue dashed line refers to a sequence with (K) = 155 and the blue 
pentagon shows the equilibrium model having as central density the time-averaged central 
density of the stable solution (i.e. model "C" in table Q. All models "A, B, C" have the 
same rest-mass to a precision of 0.4% and the arrows show how the collision corresponds to a 
number of transition in the space of configurations. 



with the ones obtained from the linear perturbation theory for the corresponding TOV model. 
This is reported in table [2j which collects the oscillation frequencies as computed from a 
Fourier analysis of the central rest-mass density of the largest subcritical solution. Because 
these frequencies with their error-bars are within the expected ones, we cannot confirm the 
claim made in JT21 that the frequencies of the critical solution are 1 or 2 orders of magnitude 
smaller than the corresponding equilibrium ones. Rather, we conclude that the metastable 
critical solution is indeed only a perturbed TOV star. 



3.3. On the critical exponent and its fluctuations 

The theory of critical phenomena predicts a precise scaling relation between the survival time 
of the nearly-critical solutions, namely the time over which a metastable equilibrium exists, 
and the distance from the critical solution. At lowest order this scaling relation is a simple 
power-law of the type 

r oq ~ -7 In \p c - p* | + const . (9) 

Following (T), we measure the survival time by considering the function £(£) = (a — 
a*) /a*, where a is the lapse function at the coordinate origin of a given simulation and a* is 
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Table 2. Eigenfrequencies of the modes of the critical solution during the metastable phase 
of the dynamics as computed from the evolution of the the central rest-mass density. Also 
indicated are the first overtones of the fundamental mode for a TOV star constructed with 
K and p c equal to the time-averages of (K) and (p c ) during the metastable phase. Despite 
the large uncertainty due to the very short integration time, the match between the two set of 
eigenfrequencies is very good. 



Mode 


Observed freq. 


(kHz) 


TOV freq. (kHz) 


HI 

H2 
H3 
H4 


4.8 ±1.2 
7.2 ±1.2 
9.6 ±1.2 
12.0 ±1.2 




3.95 
6.86 
9.42 
11.85 
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Figure 6. Left panel: Survival time of the metastable solution plotted against the logarithm 
of the difference between the initial central density of the stars and the critical one. The red 
triangles represent the data points obtained with a grid resolution of h = 0.1 and the blue 
squares represent the ones obtained with h = 0.08. The black dashed line represents the 
expected power-law scaling. Right panel: Harmonic fluctuations in the critical exponent 7 as 
obtained after subtracting the power-law scaling from the data points computed with h = 0.1 
(red triangles). Indicated with a blue dashed line is the fitting sine function. 



the lapse of the best numerical approximation of the critical solution. We set r oq to be the first 
(coordinate) time at which £( r cq) > 0.05; as discussed in (TJ, the determination of the critical 
exponent does not depend sensitively on this cut-off time. Finally, we compute 7 performing 
a linear least-square fit of (|9]i on the data points. The results of this process are shown in the 
left panel of figure [6j where we report with red triangles the values of r cq as computed from 
about 60 simulations having different initial central density. As we will discuss below, such 
a large number of data points is necessary not only to measure accurately the exponent 7, 
but also to determine whether nonlinear corrections to expression (O should be considered. 
In this way we have computed the critical exponent to be 7 = 10.004, which agrees within 
8.4%, with the value computed by JT1. As a further validation, we have computed the critical 
exponent also for a (smaller) set of simulations carried out at a higher resolution (i.e. h — 0.08 
vs h = 0.1) and these are shown as blue squares. These higher-resolution simulations predict 
a critical exponent of 10.303, thus with a difference of 2.9% from the lower-resolution ones. 
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A more careful analysis of the data for the survival time reveals that relation (0 is well 
reproduced by the data, but also that the latter show additional, fine-structure features which 
are are not accounted for in ([9j. In particular, it is apparent already at a visual inspection that 
the critical exponent also shows a periodic change as the solution approaches the critical one. 
This is highlighted in the right panel of figure [6] where we show the deviations of the data 
from the relation (|9), and where it is apparent that these deviations are essentially harmonic 
in the range in which data is available (for simplicity and to increase the precision of the fit, 
we have considered only the low-resolution data, for which more simulations are available). 
As a result, we can correct the scaling relation (0 with a simple expression of the type 

i~ C q — -7 m \Pc - P*\ + ci sin (c 2 In \p c - p*\ + c 3 ) + const . (10) 

where c\ ~ 2.0, c 2 — 1.2 and C3 ~ 0.8. Interestingly, this fine structure of the time-scaling 
relation has been observed also in in the critical collapse of scalar fields lfT3l[T4ll . but has 
never been reported before for perfect fluids, although it may be present also in the data of HI 
(cf. their figure 4). Because the scaling relation (|9]l is derived after performing a linear analysis 
near the critical solution Q, the additional oscillation captured in expression ( fTOb is a purely 
nonlinear effect which has not been yet fully explained. 

4. Conclusions 

Critical phenomena are ubiquitous in many different branches of physics and are of great 
interest in general relativity where they are associated with phase transition of families of 
solutions. With the goal of studying the occurrence of type-I critical collapse, we considered 
the head-on collision of equal mass, nonrotating NSs boosted towards each other. After fixing 
the initial velocity of the stars, we evolved numerically a great number of configurations with 
different initial central rest-mass density using the 2D general relativistic hydrodynamical 
code Whisky2D. 

Overall, the basic dynamics of the process is rather simple: As the two NSs are accelerate 
towards each other, they collide leading to a merged object which is wildly oscillating and with 
a gravitational mass which is above the maximum mass of the initial configuration. Depending 
on whether the initial central rest-mass density is larger or smaller than the critical one, the 
metastable solution either collapses to a black hole (supercritical solutions) or expands to a 
new stable stellar solution (subcritical solutions). Exploiting the accuracy of the code and the 
large set of initial configurations considered, we were able we are able to identify the critical 
value for the central density with an accuracy of 1 1 significant digits, a level of precision never 
achieved before in this type of study. 

Although overall more accurate, much of the results reported here about the critical 
behaviour coincide and confirm those found by Jin and Suen in QJE). However, we do 
differ and significantly when it comes to the interpretation of the critical solutions. More 
specifically we have shown that the head-on collision of two NSs near the critical threshold 
can be seen as series of transitions in the space of configurations from an initial stable model 
over to a metastable one with the same rest-mass but larger gravitational mass as a result of the 
conversion of the kinetic energy into internal energy via large shocks. The metastable solution 
is on the linearly unstable branch of the equilibrium configurations and thus it can exhibit 
a critical behaviour (see paper I) and either move the stable branch of equilibria (subcritical 
solutions), or will collapse to a black hole (supercritical solutions). Hence, the critical solution 
is indeed a (perturbed) TOV solution, in contrast with the conclusions drawn in lfTlfT2l 

Finally, we have also computed the critical exponent of the scaling relation of type-I 
critical phenomena and found it to agree well with the one computed by |T). However, we 
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r hoop/ 2M ADM 

Figure Al. Spacetime diagram showing the worldines of the proper stellar radius for either a 
subcritical solution R su bcrit. (black solid line) or for a supercritical solution -Raupercrit. (red 
dashed line). Both quantities are expressed in units of twice the ADM mass, so that the vertical 
dotted line at 1 represents the threshold below which the hoop conjecture would be violated 
(the time coordinate is the proper time of an observer comoving at R). Also indicated with a 
blue long-dashed line is the worldline of the apparent horizon which is relevant in the case of 
the supercritical solution. Note that the metastable solution is always outside its "hoop". 



have also found that, superimposed with the standard power-law, the scaling law shows a 
fine structure in terms of a periodic fluctuation. These wiggles in the critical exponent have 
already been observed in the case of type-II critical collapse of massless scalar fields 11311141 . 
but were never reported before in the case of perfect fluids. 
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Appendix A. Metastable solutions and the hoop conjecture 

As discussed in the main text, near the critical solution and during the metastable stage, the 
evolution of all the hydrodynamical and field variables is essentially the same (to a precision 
which is proportional to the distance from the critical solution) for both subcritical and 
supercritical solutions. It is therefore natural to ask whether the metastable solution, which 
can either produce a black hole or an expanded star, is ever compact enough to violate the 
hoop conjecture. We recall that the conjecture states that a black hole is formed if and only if 
a "mass" M of matter is confined in a hoop of radius which is in every direction is smaller or 
equal to the corresponding Schwarzschild radius. Stated differently, a black hole is expected 
to form if the matter is enclosed in a hoop of proper circumference C < AttM. 



Critical Phenomena in Neutron Stars II: Head-on Collisions 15 

The hoop conjecture, as it was originally formulated, is not meant to be a precise 
mathematical statement l24ll . even though it can be made precise under particular 
circumstances |25-27j. Most importantly, it is meant mostly as a qualitative description of 
the gravitational collapse due to the compression of matter. In particular the conjecture leaves 
much freedom in the definition of both C and M. In view of this freedom, we have taken C to 
be the proper circumference of a coordinate circle enclosing 95 % of the mass of the system 
(in order to exclude the extended low-density regions) and used the ADM mass within C as a 
local measure of the mass (This use of the ADM mass is strictly speaking incorrect as such 
a mass is well defined only at spatial infinity; however it provides an approximation which is 
adequate for the quality of the arguments made here). 

In figure lATl we show in a spacetime diagram the worldines of the proper stellar radius for 
either a subcritical solution i? S ubcrit. (black solid line) or for a supercritical solution -R S upercrit. 
(red dashed line). Both quantities are expressed in units of twice the ADM mass, so that the 
vertical dotted line at 1 represents the threshold below which the hoop conjecture would be 
violated (the time coordinate is the proper time of an observer comoving at R). Also indicated 
with a blue long-dashed line is the worldline of the apparent horizon which is relevant in the 
case of the supercritical solution. Quite clearly the metastable solution is always outside 
its "hoop" (indeed about twice as large) and when the supercritical solution crosses it is to 
produce a black hole as the conjecture predicts. It remains unclear how these worldlines 
would evolve if we had considered stars with larger boosts or that are are far from the critical 
solution. These questions will be addressed in our future work on the subject. 

As a final remark we note that if we had used the interpretation of the head-on collision 
as a transition between a stable and a metastable TOV solution {i.e. the one summarized in 
figured), then it would have been rather obvious that the hoop conjecture cannot be violated: 
a TOV star has a surface which is always outside of its Schwarzschild radius. 
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